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Abstract 

The STARWALL/CAS3D/OPTIM code package is a powerful tool to study the 
linear MHD stability of 3D, ideal equilibria in the presence of multiply-connected 
ideal and/or resistive conducting structures, and their feedback stabilization by 
external currents. Robust feedback stabilization of resistive wall modes can be 
modelled with the OPTIM code. Resistive MHD studies are possible combining 
STARWALL with the linear, resistive 2D CASTOR code as well as nonlinear MHD 
simulations combining STARWALL with the JOREK code. 

In the present paper, a detailed description of the STARWALL code is given and 
some of its applications are presented to demonstrate the methods used. Conduct¬ 
ing structures are treated in the thin wall approximation and depending on their 
complexity they are discretized by a spectral method or by triangular finite ele¬ 
ments. As an example, a configuration is considered consisting of an ideal plasma 
surrounded by a vacuum domain containing a resistive wall and bounded by an 
external wall. Ideal linear MHD modes and resistive wall modes in the presence 
of multiply-connected walls are studied. In order to treat the vertical mode self- 
consistently the STARWALL code has been completed by adding the so-called Liist- 
Martensen term generated by a constant normal displacement of the plasma. 

The appendix contains the computation of the 2D Fourier transform of singular 
inductance integrals, and the derivation of an asymptotic expansion for large Fourier 
harmonics. 


Keywords: resistive wall modes, vertical modes, energy principle, multiply-connected 
wall structures, variational method, finite element method 
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1 Introduction 


The STARWALL code has originally been developed for the numerical treatment of re¬ 
sistive wall modes (RWMs). It has been applied to many different physics problems via a 
coupling to linear and non-linear MHD codes already (see also the outlook in Section 6). 
The present paper concentrates on STARWALL itself, in particular on the mathematical 
and numerical methods used. 

After briefly introducing resistive wall modes, this introductory part explains the main 
features of the STARWALL code and its interaction with other codes. Feedback stabiliza¬ 
tion studies of RWMs are described, and the recently added Liist-Martensen term which 
is important for axisymmetric instabilities is introduced. Finally the outline of the rest 
of the paper is given. 

External kink unstable tokamak equilibria which are fully stabilized by an ideal conducting 
wall sufficiently close to the plasma remain unstable in the presence of a a realistic wall 
because of its hnite resistivity. These instabilities. Resistive Wall Modes (RWMs), grow on 
the time scale of the magnetic field diffusion through the resistive wall. With the growth 
rate of the RWM being typically orders of magnitude smaller than that of the kink mode 
in the absence of a wall, stabilizing the modes with an active feedback system becomes 
feasible. The topical review on stabilization of the external kink and the resistive wall 
mode by Chu and Okabayashi [1] gives a comprehensive overview of the existing literature 
on this topic. 

Beside STARWALL - presented here - several other codes have been developed and used 
to study RWMs: VALEN [IE], DCON coupled to VACUUM m, MARS-F |S|, CarMa [Sj. 
In the VALEN code the plasma state is approximated by a single unstable eigenfunction, 
whereas in MARS-F and CarMa the stability of 2D equilibria in the presence of 3D 
structures is treated. STARWALL is the only code which can be applied to 3D equilibria 
with general 3D wall configurations. Both, VALEN and STARWALL use a thin wall 
approximation. 

The STARWALL code is part of a comprehensive code package: 3D equilibrium NEMEC 
code HlEli, coordinate transformation COTRANS code, 3D ideal MHD stability CAS3D 
code [lO], vacuum and RWM stability STARWALL code, and the feedback optimization 
OPTIM code [H]. 

An adapted version of the STARWALL code has also been coupled [12] to the non-linear 
MHD code JOREK [T3] to include resistive wall effects in the simulations and has already 
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been applied to a variety of different physics problems. The CASTOR3D code is currently 
under development la, for linear stability studies of 3D equilibria including 3D resistive 
wall effects described in the STARWALL formalism. Some additional information on 
JOREK-STARWALL and CASTOR3D is provided in Section 6. 

Computations of external modes with ideal conducting wall conhgurations are possible 
with the CAS3D code including the perturbed kinetic energy where the vacuum energy 
term is provided by the STARWALL code. In case of slowly growing resistive wall modes, 
the plasma inertia can be neglected such that a problem of hrst order in time is ob¬ 
tained |1]. These modes can be studied with the STARWALL code. Also, the feedback 
stabilization of RWMs can be investigated as described in the following. The feedback 
procedure can be divided in two parts, the open-loop and the closed-loop problem. In 
the open-loop part, a complete set of eigenfunctions of the plasma-resistive-wall system is 
determined. The feedback coils can be included in the resistive wall conhguration passive 
resistive elements without external voltage applied. 

The closed-loop part consists of feedback logics which calculate the optimal voltages to 
be applied at the feedback coils based on signals which are produced by sensor coils 
at appropriate positions. The feedback code OPTIM m implemented for this purpose 
achieves robust control by considering all unstable, but also a set of stable modes to assure 
that modes are not driven unstable by the active feedback. 

RWM kink modes and their feedback stabilization have been studied for ITER and AS- 
DEX Upgrade tokamak-type conhgurations with realistic wall structures and for 3D quasi- 
axisymmetric equilibria [la usi El na [19]. 

To study external axisymmetric modes (toroidal harmonic n = 0) the STARWALL code 
has been completed by implementing the so-called Liist-Martensen term [20]. The vac¬ 
uum contribution of an external mode is determined by the normal displacement at the 
plasma boundary. For ^ const the vacuum contribution is the solution of a Neumann- 
type problem with prescribed normal component B„. The Liist-Martensen term is the 
plasma perturbation for = const where the perturbed is tangential at the plasma 
boundary. In that case a net-toroidal and net-poloidal magnetic hux is induced in the 
vacuum region bounded by the plasma boundary and an external ideal conducting wall or 
at least by ideal toroidal held coils. That is, for vertical mode studies the identical phys¬ 
ical conhguration should be used by which the free-boundary equilibrium was generated. 
Results of Vertical Displacement Events (VDEs) taking into account the Liist-Martensen 
term and their coupling to higher toroidal harmonics via a three-dimensional resistive 
wall are presented for an AUG-type conhguration in Section 5. 
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The rest of this paper is organized as follows. In Section 2, the physical problem to be 
solved for vacnnm, resistive wall, and ideal plasma is dehned. In Section 3 the variational 
procednre to solve this problem nsing a spectral discretization is described, and in Section 
4 the method nsed for triangnlar hnite elements is introdnced. Section 5 contains an 
example for applying STARWALL to realistic geometries: Linear stability is investigated 
for an ASDEX Upgrade-type conhgnration inclnding a 3D resistive wall with holes. In 
Section 6, a brief ontlook to fntnre applications of the STARWALL code for linear and 
non-linear MHD problems is given. 

In Appendix A the variational method is explained comprehensively, and Appendix B 
contains the dehnition of the indnctance matrices. In Appendix C the snbtraction method 
for the treatment of the singnlar matrix elements is described. Appendix D contains the 
compntation of the 2D Fonrier transform of singnlar indnctance terms already pnblished 
in [21] and [22] . The treatment of nnstable recnrrence relations by a method rnnning the 
recnrsion in the backward direction has been added. Fnrthermore, the derivation of an 
asymptotic expansion for large valnes of the 2D Fonrier harmonics m, n for the singnlar 
Fonrier integrals is given. 

2 The vacuum contribution 


The STARWALL procedure will be explained considering an ideal plasma equilibrium 
in the presence of a multiply-connected resistive wall and an external ideal conducting 
(superconducting) wall (see Fig. 1 and Fig. 2). For the resistive wall the thin-wall ap¬ 
proximation is used [23]. Assuming the RWM to be sufficiently slow, the kinetic energy 
of the plasma perturbation can be neglected. In analogy to the ideal energy principle [21] 
a variational technique can be applied. The Lagrangian 

C = W/tS) + l f df (1) 

V Sj. 

consists of the potential energy of the plasma perturbation and the contribution of the 
vacuum domain given by a surface integral at the plasma-vacuum interface {Si = plasma- 
vacuum interface, ^ = displacement vector, Bo= equilibrium magnetic held, B= perturbed 
vacuum magnetic held, n = exterior normal on Si). The perturbed magnetic held B is 
uniquely determined by the normal component n ■ ^ of the displacement ^ at Si. 
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The potential energy of the plasma perturbation pa - provided by the CAS3D code - is 
given by 

^ ll + rp(v■■ nfl, (2) 

V s± 

c = V X X Bo) + (^ ■ n)j X n, A = 2(j x n) ■ (Bq ■ V)n, 

Bo = Vs X V(Fp u — m), 

where s, u, v are flux coordinates: s : 0 < s < 1 is the normalized toroidal flux, and 
{u,v) : 0 < u,v < 1 are poloidal and toroidal magnetic coordinates on the surfaces. 
Fp, F!p are the derivatives of the poloidal and toroidal flux with respect to s. The Fourier 
expansion of the displacement vector reads 

mf,nf 

^{s,u,v) = E ^^(s)mn sin 27r(mn + nv) + ^cos 27i{mu + nv). (3) 

m=0,n=0 
m>0,n = —n/ 

With respect to the flux coordinate s, the Fourier harmonics ^s{^)mnAci^)mn are dis¬ 
cretized using a hnite element method. 

The perturbed magnetic held B has to satisfy Maxwell’s equations 


B = VxA, Vx(VxA)=0, V-A = 0 


(4) 


with boundary conditions for the vector potential A at the plasma-vacuum interface, the 
resistive wall and the external superconducting wall. 

On the resistive wall the boundary condition follows from Faraday’s and Ohm’s law: 
E ^ = 0, crE = J . Assuming a time dependence e'*'*, one gets in the thin wall 
approximation the boundary conditions for the perturbed vector potential 

{ -(n-^)Bo : on Si (plasma-vacuum interface) 

-^nxj 2 : on S '2 (resistive wall) ( 5 ) 

0 : on S '3 (external conducting wall) 


where j 2 is the current in the resistive wall, n is the exterior normal unit vector on the 
surfaces, ad is the surface resistance of the wall {a = conductivity, d = wall thickness) 
and Bq is the equilibrium magnetic held. 


With the contravariant normal component Vs ■ ^ used in the CAS3D code the boundary 
condition at the plasma-vacuum interface reads 


Vs X A = —(Vs ■ ^)Bo, 


Vs 


(r^ X r„) 


n, 


Fp Fp Yu 

(r^ X rj ■ 


( 6 ) 
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with the equilibrium field Bq in magnetic coordinates and r„ := r„ := := 

Multiplying ([^ with n x r^j and n x Fv one gets 


ri,« - A 

= Ff Vs-t 

- A = 

-F'p Vs-^ 

r2,u - A 

1 

r2,v - A = 

1 „ : 

adj 

7 J 2 

ad^ 

r3,« ■ A 

= 0, 

rs,.; ■ A = 

0 


on Si 

on S 2 (7) 

on S's 


The boundary conditions for the perturbed magnetic field B are obtained by taking the 
divergence: [V ■ (Vs x A) = A ■ (V x Vs) — Vs ■ B], At the plasma-vacuum interface Si 
one gets 

Vs-B = Bo-V(Vs-^) ( 8 ) 

or 

{ F^ (Vs ■ + Fp (Vs ■ ^)u ■ on Si (plasma-vacuum interface) 

^{r 2 ,v ■ h,u - r 2 ,u ■ h,v) ■ on S '2 (resistive wall) (9) 

0 : on S '3 (external conducting wall) 

with N = X r„ to be taken at the surfaces. 

Prescribing the normal component Vs-B of B one has to solve a Neumann-type problem 
except for Vs - ^ = constant where the normal component of B vanishes. To get the 
boundary condition for this case one has to go back to the boundary condition for the 
vector potential A. 

The solution for the vector potential A can be generated by surface currents on the 

plasma-vacuum interface Si, the resistive wall S '2 and the external ideal conducting wall 

S 3 



with divergence-free surface currents derived from current potentials 


ji = iij X V<h*, rij = 


r, ^ = 1,2,3 


where rij, z = 1,2,3 are the exterior unit normal vectors. 

For toroidally and poloidally closed tori the current potential is given by 


( 11 ) 


= if u + if V + (pfu, v), z = 1, 2, 3 


( 12 ) 
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where If, If are net-toroidal and net-poloidal currents on the torus. They play a role 
only if the diplacenient ^ contains a component Vs ■ ^ = constant [20]. The (pfu^v) are 
single-valued potentials which can be expanded in Fourier space for smooth tori 

mc,nc 

v sin27r(mn-|-?w)-(-0^(m,n) cos27r(mM-|-?7,n), i = l,2,3 (13) 

m= 0 ,n = l 
Tn>0,Ti= —ric 


The surface currents have to be determined such that the boundary conditions for A on 
Si, i=l,2,3 are satished. That will be carried out by means of a variational procedure. 

3 Variational method 


To treat cases with multiply-connected resistive wall configurations, a finite element 
method has been applied using a variational procedure [23|. One introduces the La- 
grangian 




i=l k=l 


Ji ■ Jfc 

r* - rfc 


S2 


Sk 


Si 


(14) 


In appendix A it is shown that the hrst variation SCy = 0 gives the boundary conditions 
(A.10),(A.11),(A.12). Assuming a smooth plasma-vacuum interface being identical with 
the plasma boundary the current potential and the normal component of the displacement 
vector can be expanded in Fourier space 


mf,nf 

^ ■ Vs = .^s(m, n) sin 27r(mM-|-nn)-f |c(’^,''^) COS 27r(m-u-|-?w) (15) 

m = 0,n=0 
7n>0,n = — n ji 


so that one gets for the last term of the Lagrangian Cy 


y#! (nr^)ni-(jixBo) = -{Fflf + Fflf) i{0,0) (16) 

Si 

mf,nf 

-I- 'y^7i{nFf -|- mFp)(,^^(m, n) — 

m=0,n = l 
m>0,n = —ny 

being linear in and determining the boundary condition at the plasma-vacuum in¬ 
terface. The variables of the Lagrangian Cy are the net-currents If ,lf and the Fourier 
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harmonics of the single-valued potentials = 1,2,3. Replacing in the currents 
using 0 one gets 


3 3 

+ ^ ^2 N 22 <h 2 + $7 


i=l k=l 


27 


(17) 


with 


*^7^ = If, 01(0,1), «/J, 0i (0,1),..., 0-(m/^, u/J] 

= [ 6 ( 0 , l),...,|.(m/,n/),|c( 0 , 0 ),...,^c("^/,?^/)] (18) 

and T = transpose. 

Varying Cy with respect to <l)j, z = 1, 2, 3 one gets a set of linear equations for . 

Mpp <l>p -|- Mp^ -|- Mpt ^ 

M.UJP <l>p + -t- = 0, (19) 

Mtp <l>p -|- <!>.„, -|- Mtt = —Nit “hi 

7 

(note: the index notation has been changed p = 1 - plasma-vacuum interface, t = 2 - 
resistive wall, w = 3 - external conducting wall) 

The vacuum contribution in ([T| is given by 

Ws = ^ f du du(^ ■ Vs)(Fp -f F0 ri,^) ■ B (20) 

V s± 

with the magnetic held generated by the currents on the surfaces Si, i = 1,2,3 


3 r 


Bfri = 


-V 

dvr 


2=1 




du' dv 


>Si 


— ^'13 


r — r 


-HV / (df' -V' 


'Si 


r — r 


v) 0 ( 


Inserting (21) into (20), Wg splits into two terms 

1 ^ 


m'. = 2E(V‘’ + W''; 


( 2 )- 


( 21 ) 


( 22 ) 


2=1 


The contributions for z = 1 are given by 


WA ^ ^ y '^^(^ ■ '^®)( J i.F'p + Ft rpi;) x {if - if riy) 


(ri - 6 ) 
|ri - r;|3 


27i{F'p If + F!p If 


ivf = j dudv[{F'p^ + F!r^){i-Ws)\{2T: 4>{u,v) 

J Si 


(j- _ J-M 

-f 4 dudv'{viy X ri,„/) ■ -6 - f- 0 i(zz',z;') 


Si 


(23) 












For i = 2,3 the contributions to (22) are given by 




( 1 ) _ 


= — I du dv{i ■ Vs) / du'dv'{F'p ri,^) x {if - if Yiy) ■ 


(ri - r •) 


'Si 


's. 


— r'|3 


iy/2) ^ 1 

Att 


du dv 


'Si 




'Si 


du'dv'{ri^y> X Tiy) ■ ^ (j)[ (24) 




The u, n-integrations have to be performed infinitesimal exterior of the plasma-vacuum in¬ 
terface Si and the singular integrals wf \ wf'^ are treated using a subtraction method(see 
appendix B). 


With the dehnitions (18) the vacuum contribution Wg reads 

n',= i(r 4.,+•*„+r M^, •*,) 


(25) 


The elements of the matrices and are given in appendix B. 

The variables of the system considered are the Fourier harmonics of the displacement 
vector the current potentials on the plasma-vacuum interface, the resistive wall and 


the external wall. The set of equations (19) determines the current potentials for given 
normal component ^ of the displacement vector on the plasma-vacuum interface. 

The system is closed by the equation derived from the Lagrangian ([^. Varying the 
discretized functional ([^ one obtains 


( ^ ) ( dp -t- dt 

The variables are denoted as follows: X = ^ ■ ■ ■} consists of all 

components of the displacement vector except for the Fourier harmonics of the normal 
component at the plasma boundary which are denoted by ^ = {Vs ■ ^ 5 (l)mn, • • •, Vs ■ 

Solving the set of equations for ^ one gets 



V ^ix J 


with 


^ casSd /V, /V, /V, 

^ (Mj, % + + M|, 4,) 


casSd 


(27) 

(28) 
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Eliminating ^ in (19) one obtains 

~ M 


with 


Mpp 


M 


pw -|- IVIpj — 0 , 

<l>t = 0 , 

Mu < 1 , = - 


Wp ^ p 

Mtp $„ + 


ad'j 


Ntt 


casSd 


Mpp = (Mpp + M. M-M,J 


pi a ip' 

casSd 

-—- cas3d 

Mpi = (Mpt + 


Mp^ = (M, 


pw 


(29) 


(30) 


Finally one gets a set of linear equations defining an initial value problem of first order in 
time. Assuming a time dependence the normal modes of the system are obtained by 
solving the generalized eigenvalue problem 
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M 


tt 


Mtp M: 


tw 


'( 


Mpp M, 


pw 




-1 


M 

M 


pt 


wt 


= —-Ntt <ht (31) 
ad 


4 The finite element method 


For the finite element procedure the conducting wall is discretized into triangles. The 
position vector of a triangle is given by 

r = n+ « r 2 ,i+/d rsj, q(> 0, /3>0, 0<a + /?<l, ri^k = ri-rk, z, A: = 1, 2, 3 (32) 
with the vertices numbered anti-clockwise (ri,r 2 ,r 3 ) and r 4 = ri. 

Assuming the surface-current density to be constant on the triangle the current can be 


written as 


JA = 


01 1'2,3 + 02 fS,! + 03 I‘p2 




or Ja ^ ^ 0j Oi) Oi „^ikl 

^ 2 |r 2 ,i X r 3 , 2 | 


T 2.1 X r3,2 . , 

where the 0* are the values of the current potential at the vertices of the triangles. 
The contribution of a pair of triangles to the functional C is given by 

1 


(33) 


Ca'a =}a' ■ Ja^ [ df [ df — 

Stt Ja' Ja I r' 


(34) 
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^A'A ^ ^ ^ 4^iLik(l^k'! 


i,/c=l 


fj^'fjiAr 


The fourfold integral Li^k can be computed analytically. One gets a sufficiently good 
approximation by performing two integrations analytically and the remaining two inte¬ 
grations numerically: 

One obtains for the twofold integral 

df \__ \_,i = ^|r^+i,i|(a^ In (35) 


r — r 


i=l 

hi ( arctan 


C + C-l 

tti df 


if + hf + If hi 


arctan 


a,- d. 


* 


+ hf + li hi 


with 


It = 


ki+i - r| _ |ri - r| ^ _ (r^+i - r) ■ r^+i^i _ (r^ - r) ■ Yi+i^i 

km.l ’ * “ km.l ’ * “ km.P ’ * “ km.P 


(Fi - r) X ■ n (ri-rj-n r2,i X 

di = -:- 77. -, hi = -;-;-, 11 = 


ki+i,il 


ki+i,il 


T2,1 X 13 ,i| 


The remaining two integrations are done numerically using &Ng = 3-point or Ng = 7-point 
Gauss quadrature formula for a triangle domain 


df / df 


r — r 




dfi 


2=1 


ri + r2,i 0 + rs.i Vi “ 


(36) 


with (i, rji, and Wi listed in Table I. 


i 

0 

Vi 

Wi 

1 

1/3 

1/3 

270/2400 

2 

(6+ v^)/21 

(6+ yT5)/21 

(155 + yi5)/2400 

3 

(9-2yi5)/21 

(6+ yT5)/21 

(155 + yi5)/2400 

4 

(6+ v^)/21 

(9 - 2v^)/21 

(155 + yi5)/2400 

5 

(6- yi5)/21 

(6- v^)/21 

(155 - yT5)/2400 

6 

(9 + 2v^)/21 

(6- v^)/21 

(155 - v^)/2400 

7 

(6- yi5)/21 

(9 + 2yT5)/21 

(155 - yT5)/2400 


i 

0 

Vi 

Wi 

1 

1/6 

1/6 

1/6 

2 

2/3 

1/6 

1/6 

3 

1/6 

2/3 

1/6 


(a) (b) 

Table I Gauss quadrature with weights Wi at points Q and rji 
for a. Ng = 7 point (a) and Ng = 3 point (b) formula 
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For the self-inductance of a triangle one gets 


-hjfc Gj ■ G/j 

with L = |r 2 ,i| |r3,2| + kgsl- 


|r2,i X rg^i 
127r 


2 3 


E 


j=i 


ln( 


L 2 |rj+ij| 


(37) 


The magnetic held Ba produced by a constant current j a on the triangle can be computed 
analytically 

BA = -;^jAxv/ (38) 

47r Ja' |r - r'| 

with 


df, 


r — r 


E 

i=l 


n X + 1- +1 


ki+i,. 


In 


+ c 


(39) 


r, - r ■ n 


Ht 


(r* - r) ■ n| V 


arctan ■ 


tti di 


+ 


if + hf + ItK 


— arctan 


tti d~: 


+ dj + li hi 


so that the contribution of a triangle to the vacuum energy term (20) is given by 


Wa = 1; [ du dviVs ■ ^)(F^ri_„ Fpri,„) ■ Ba 
^ Jsi 


(40) 


The surface-current potential on a poloidally and toroidally closed surface consists of two 
multivalued secular terms determining the net-poloidal and net-toroidal current and of a 
single-valued periodic term (j). The current potential is approximated by its value at the 
numbered global nodes of the triangulated domain. There are three local nodes at the 
vertices of each triangle. The nodes are numbered anti-clockwise. One dehnes the vector 
of the nodal values at the vertices of all triangles by 


^ind 

$(*) = + * = (41) 

i=i 


where Nt is the number of triangles and Nind is the number of independent variables 0(j) 
at the global nodes. A global node belongs to several triangles. U{i) and V{i) are the 
nodal values of a unit net-toroidal and net-poloidal current distribution which can be 
arbitrarely chosen. The so-called connectivity matrix H relates the global nodes to the 
local nodes of the triangles. If there are holes in the wall the current potential has to 
satisfy the boundary condition <F = const, along the edges: the nodes at the vertices of 
the triangles along the edges of a hole have the same value reducing the number of the 
independent variables. That has to be incorporated in the connectivity matrix. 
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The matrix elements of a Lagrangian £5 = are 


M,,_ 


Sz—3+/c, 3z'—3+fc' — 


% ■ 


An 


df / df 


Til - Fi 


fc = 1,2,3, i = ,Nt 
k' = 1,2,3, z'= !,■■■ ,Nt 


(42) 


The contribution of a triangulated (resistive) wall to the Lagrangian is then given by 
1 


Cs = Mtt $, Mu = G^MG, G = {U,V,H), $ = | I 


p 


(43) 


where <h is the vector of the current potential of the independent variables. The con¬ 


tribution to the vacuum energy matrix of the wall is obtained with (38),(39),(42) and 


(43). 


= M^^G, (44) 

where the matrix elements of are the contributions to the Fourier harmonics of 

Vs ■ ^ for the nodal values of the vertices of all triangles. 


One possible choice for the net-current contribution for a closed torus without holes can 
be dehned as follows: On a rectangular mesh with poloidal and n„ toroidal meshpoints 


'^k) 


/i — 1 fc — 1\ 

V ’ n„ J 


i = 1 , • • • ,nu + l 

’ A: = 1 , • • • ,n^ -f 1 


(45) 


the positions of the 2nuny triangles are given by 


= r{ui ,Vk+i), ri^i = r{ui+i,Vk ), (46) 

Z — 1 

r;_i,2 = r{ui+i,Vk+i), ri^2 = r{ui ,Vk ), I = 2{i + nu{k-l)), ’ 

k = l,nv 

r;_i,3 = r{ui ,Vk ), fi^3 = r{ui+i,Vk+i) 


where I labels the triangles. There are (n„ -f- l)(n^ -|- 1) nodes for the secular terms being 
proportional to and . The values of these current potential terms at the nodes can 
be chosen as 


U{i) = Ui 

V{j) = Vk 


i = i + {nu + 1 ) {k - 1 ), 


i = !,■■■ ,nu + l 
k = 1,--- ,n„ + l 


(47) 


while there are nodes for the periodic current potential 0. 


A resistive wall closed only poloidally or toroidally is topologically a cylinder. In that case 
the current potential becomes single-valued. For a poloidally (toroidally) closed wall the 
net-toroidal (net-poloidal) current vanishes and the net-poloidal (net-toroidal) current is 
given by the difference of the current potential value at the two boundaries of the wall. 
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5 Application 


Numerical results are presented for an ASDEX Upgrade-type test equilibrium which is 
unstable with respect to external modes. The plasma equilibrium properties are: major 
radius Rq = 1.64 m, plasma current Ip = 0.98 MA, monotone g-pro£le with qaxis = 1-46 
and qboundary = 5.26, vacuum magnetic held strength Bq{Rq) = 2.43 T, and volume 
average beta < /3 >= 2.58%. A poloidal cross-section of the hux surfaces, as well as the 
pressure and g—profiles are shown in Figs la-c. In Fig. la the positions of the plasma 
boundary (red), the resistive wall (blue), and the toroidal held coils (brown) are sketched. 
At the low held side the plasma boundary extends to i? ~ 2.14 m, while the resistive wall 
is localized at i? ~ 2.23 m, so that the plasma-wall distance amounts to AR ^ 9 cm at 
this position. 




0.0 0.2 0.4 0.6 0.8 

Ptor 


1.0 


Fig. 1 : (a) Poloidal cross-section of Hux surfaces, plasma boundary (Si), resistive wail 
(S 2 ), and toroidal Held coil (S 3 ), (b) pressure proHie, and (c) q-proHie of the ASDEX 
Upgrade-type test equilibrium (ptor = \fs)- 
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Figure 2a shows a 3D view of the preliminary wall design of ASDEX Upgrade [26]. In 
Fig. 2b the whole conhguration composed of plasma, resistive wall, and 16 toroidal held 
coils is presented. The latter are modeled by inhnitely thin bands of hnite width. 



-\- \ -1-1-!-1 

• 3 - 2-10 1 2 3 


(a) (b) 

Fig. 2: (a) Preliminary design of the resistive wall, and (b) plasma with resistive wail 
and toroidal Geld coils. 

Without wall the equilibrium is unstable with respect to u = 0, 1 and 2 modes (n > 2 are 
not considered). However, it can be stabilized with an ideal conducting wall sufficiently 
close to the plasma (see Fig. la) 

In case of a hnite wall conductivity the plasma is unstable on a resistive time scale. A 
surface conductivity ad = 2.8 ■ 10^ S was used with a being the specihc conductivity and 
d the thickness of the wall. 

In Fig. 3a the eigenvalue 7 of the n = 0 mode is plotted versus the plasma-wall distance 
^shift- The eigenvalue is computed with the CAS3D code including the perturbed kinetic 
energy term. The m = 0,n = 0 harmonic causes a current in the toroidal held coils 
producing a net-toroidal held in the region between the plasma boundary and the coils. 
In Fig. 4 the m-harmonics of the n = 0 eigenfunction are plotted for the no-wall limit 
case (7 = 510481 1/s). There, the m = 0 harmonic only makes a very small contribution. 


15 






































Fig. 3 : Growth rate of the n=0 mode in dependence of the wall position. The latter is 
presented by the radial shift of the wall with rshift = 0 being the planned position of the 
wall in ASDEX Upgrade. The hatched area marks the plasma with its boundary lying 
approximately 9 cm inside the unshifted wall. The plots show the growth rates assuming 
an ideal, and a resisitve wall, respectively. The plot on the right-hand side presents an 
enlargement of the left-hand side plot, showing the resistive wall growth rates and the 
stabilizing ideal wall position at rshift = 0.16 m in detail. 



I/Os 


2/0 s 


n.' 0 c 


3/0 c 

2 / 0 . 


Fig. 4: Eigenfunction of the n = 0 mode represented by the cosine and sine Eourier 
harmonics, ^f the normal component of the displacement vector ^ as function 

of s. Here and in the following plots, the largest harmonics are marked by their poloidal 
and toroidal indices, m/n. The s and c attached to these numbers characterize the sine 
and cosine harmonics, respectively. In this plot, additionally, the m=0, n=0 harmonic is 
marked by the orange line. 
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n=0,l,2, 7 = 24 1/s 


n=0, 7 = 26 1/s 



n=l, 7 = 50 1/s 



n=2, 7 = 314 1/s 



n=2, 7 = 276 1/s 



(a) 



n=0,l,2, 7 = 51 1/s 



n=0,l,2, 7 = 386 1/s 



(b) 


Fig. 5 : Resistive wall surface currents belonging to the 5 unstable {n = 0,1,2) multi¬ 
mode eigenvalues (b) and to the unstable single n = 0, n = l, n = 2 eigenvalues (a) are 
shown. The current pattern of the multi-mode cases can be uniqely related to the single 
n solutions. 


The axisymmetry is broken by the multiply-connected ASDEX Upgrade resistive wall 
so that, in principle, all n-harmonics contribute to an eigenmode. In Fig. 5b surface- 
current lines of the induced wall currents are shown for the 5 unstable eigenmodes where 
n = 0,1,2 toroidal harmonics have been taken into account. For comparison, in Fig. 5a 
surface current are shown obtained from single n = 0 , n = 1 and n = 2 toroidal mode 
computation. Comparing the eigenvalues - also quoted in Figs 5a-b - one can uniquely 
relate the single n modes to the n = 0,1,2 eigenmodes. That is, each eigenmode is 
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dominated by one n-harmonic. The asymmetry of the wall geometry leads to a significant 
splitting of the eigenvalues and coupling of the toroidal harmonics. In Figs 6a-b the two 
orthogonal eigenfunctions for the [n = 2) case with different eigenvalues (see Fig. 5a row 
4 and 5) are shown. In order to get a sufficiently good resolution in the neighbourhood 
of rational surfaces a non-equidistant radial grid has been implemented in the CAS3D 
code. Figure 7a shows the coupling of the toroidal n = 1, n = 2 harmonics. In Fig. 7b the 
enlargement of the 3/2-harmonic region demonstrates the improved resolution obtained 
by the accumulated mesh at the g = 3/2 surface. 




(a) (b) 

Fig. 6: Eigenfunctions of the n = 2 mode in presence of a resistive wall and ideal toroidal 
field coils. The two plots show the eigenfunction spectra for the two eigenvalues which 
correspond to the two existing orthogonal solutions. 




(a) 


(b) 


Fig. 7: (a) Eigenfunctions of the coupled n = 1, n = 2 mode, (b) Enlargement of the 
3/2-harmonic showing the radial grid accumulation around the rational g = 3/2 surface. 
To obtain the same grid rehnement with an equidistant grid 40000 grid points would he 
necessary compared to 931 non-equidistant grid points. 
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6 Outlook 


The STARWALL code is a very flexible numerical tool which can be combined with linear 
and nonlinear stability codes. The non-linear MHD code JOREK and the linear MHD 
code CASTOR3D are currently under development using STARWALL for the resistive 
wall part. 

The JOREK-STARWALL code has already been applied to nonlinear studies of resistive 
wall modes 123, vertical displacement events [28], disruptions triggered by massive gas 
injection [29] and QH-Mode ED]. Besides a current optimization and parallelization of 
STARWALL for resolving larger problem sizes, an extension of JOREK-STARWALL to 
include Halo currents is currently under development |28j . 

Furthermore, the STARWALL code is used in combination with the CASTOR_3DW code 
for stability studies of axisymmetric equilibria in presence of ideal wall structures, includ¬ 
ing physical effects such as plasma resistivity, viscosity, poloidal and toroidal rotation [31] . 
Currently, CASTOR3D is under development [H], a hybrid of the CASTOR_3DW and 
the STARWALL codes. Solving an extended eigenvalue problem (MHD equations of CAS- 
TOR_3DW and vacuum part of STARWALL), this code takes plasma inertia and resistive 
walls simultaneously into account. Besides the straight held line coordinates used for the 
plasma part so far, general flux coordinates are implemented such that linear studies of 
resistive and rotating 3D equilibria in presence of 3D resistive walls will be possible. 
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7 Appendix 


A Lagrangians 


The variational method will be explained by starting with a simple problem. For a current 
j on an ideal conducting torus one can dehne the following variational problem: determine 
the current distribution by minimizing the energy of the magnetic held produced. The 
Lagrangian - being the magnetic energy - is given by 


Cs 




(A.l) 


with the divergence-free surface current 


j = nxV<F, ^ V + u + v) 


(A.2) 


with angle-like surface coordinates u,v : 0 < u < l, 0 <n< 1 . 

Given a net-poloidal current and/or a net-toroidal current one gets the ’’Euler 
equation” by varying £5 with respect to the single-valued potential (f){u,v). 

S S S 


and with 




ru-§^6(j) - r,£6(j) 

I ft. X I 


(A.4) 


f d d f d d 

5Cs = / dti dn(r„—(50 - r^—(50) ■ A = - / dw (in50(r„—- r^ —) ■ A 

s s 

= - [du dv5(l) [r„(r^ ■ V) - r^(r„ ■ V)] ■ A 


(A.5) 


J du dvdcj) {vy X Vu) ■ {V X A) = — J du dvScj) (r„ x r^,) ■ B 
s s 


From 6Cs = 0 it follows the so-called natural boundary condition n ■ B = 0. 

There are two independent solutions. Given a net-poloidal current the magnetic held 
B vanishes in the exterior region of the toroidal surface and given a net-toroidal current 
the held B vanishes in the interior region. For the problem considered one got the natural 
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boundary condition n ■ B = 0. In order to obtain more general boundary conditions one 
has to add appropriate terms. One gets the Lagrangian Ly for the conhguration treated 


in section 3 by terms of type (A.l) and adding a term linear in j appropriate to satisfy the 


boundary condition at the plasma-vacuum interface Si and a term to fulfil the boundary 
condition at the resistive wall 82 - 


1 


3 3 




Stt 


i=l k=l 


Sk 


Vi - Ffc I 27 


J2 ■ J2 

ad 


+ / dfi (ni-^) ni-(ji X Bo) (A.6) 


S 2 


Si 


To study feedback stabilisation one has to extend the Lagrangian by adding terms for 
sensor and feedback coils and for feedback voltages given by 


ricUc 


“ to 




E I ‘if‘ I * i E I « E ] « 


Cl Cl, 


Cl 


2=1 


J» • V; 

r* - rz 


r(A.7) 


1 


TLq n -| TLq 

1=1 ^ 1=1 


feedb 

c,l 


with ric = number of coils, = coil currents, Ri =coil resistance and feedback voltages 

jj feedb 


Considering the variation of the Lagrangian Cy (A.6) with respect to <hi, i = 1, 2, 3 gives 


6 Cy 


with 


(^5lf J du dv ri^u ■ A — 5lJ' J du dv ■ A — J du dv 5(l)i (rj_^ x rj_„) ■ B j 

2=1 rt 00 


d d 


{F!^ + F'p 5Ii) j du dvVs ■ ^ +j du dv 6 (j)i {Ffp— + F'p—)Vs ■ (A.8) 

5i Si 

SI 2 fdu dv — fdu dv 5 (j) 2 - 

S 2 S 2 S 2 


6 I 2 / du dv 
^ ' 'yad 


dv ' ^du' 

^2,v ' J2,u ^2,22 ' J2,2; 
'jad 



(A.9) 
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From 5Cv = 0 one gets the boundary conditions at the plasma-vacuum interface Si 


J du dv ■ A(ri) = J du dn(Vs ■ ^), 

Si Si 

/ du dv ■ A(ri) = —Fp / du dv{'\/s ■ 


Si 


Si 


(ri,, xriJ-B(ri) = [F/,—+ F'p—){Vs ■ ^) 


at the resistive wall S 2 


du dv r 2 ^u ■ A(r 2 ) = du dv 


S 2 


S 2 


du dv r 2 ,„ • A(r 2 ) = du dv 


S 2 

{r2,v X r2,„) ■ B(r2) 
and at the external ideal wall S 3 


S 2 


r2,M ■ J 2 
'yad 

T^ 2 ,v-h 

■yad 


^2,v ' J2,u ^2,u '32,1 
■yad 


(A.IO) 


(A.ll) 


J du dv ■ A(r 3 ) = 0, J du dv r 3 ,^ ■ A(r 3 ) = 0, (r 3 ,„ x r 3 ,„) ■ B(r 3 ) = 0 (A.12) 

S 3 S 3 

The boundary conditions (A.10,A.11,A.12) obtained for A are integrals. With an appro¬ 
priate gauge A— > A -|- VT one can get the local conditions ([^. 


As shown before, on a toroidal surface there exists for an arbitrary prescribed net- 
poloidal(net-toroidal) current (/^) a so-called ’equilibrium’ current distribution 
( jjg ) generating a magnetic field with vanishing normal component on the surface. The 
held produced by such a poloidal current on Si and a toroidal current on S 3 is zero in the 
domain bounded by the boundaries Si and S 3 . With these equilibrium currents one can 
compensate the and /J contributions to the solution. Therefore one can omit these 
terms. 


B Inductance matrices 


The elements of the matrices Mjj, i,j = 

1 , 2,3 

in (17) 



/ 9^1 

9\2 

9ic' 

9is 

M — 

921 

922 

ij 

92c' 

ij 

-92 s 

iVljj 

9ci 

9c2 

9cc' 

ij 

-9js 


\ -9:{ 

-9li 

9sc' 

9ss' 
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(B.l) 








are defined by 


ij n 1 r\ ^ f J f ^ f (1 R-yx') /,/ /N 

g^ci\k,k) = — J ax J ax cos{kx) -^- j -^— cos{k x ) 

Si Sj ^ 

= “ 7 “ [ dx [ dx'^^ — ^^-—^^^j^cos{k'x') 

dTT y J I ri - r. I 

Si Sj 

dic'i^') ~ ^ ^ ^ cos{k'x') 


(B.2) 


g\\ = — dx dx 


/ ^i,v ■ ^j,v' ij _ 1 


47r 


Si Si 


r* - r . 


921 = 


I 11%,u ■ llj,v' ij _ 


Ti - r i 


47r 

1 

47r 


dx / 

'dx',''*’" 

J 

1 r* - r^- 1 

^3 


dx / 

dx' 


r,; - Fi 


with the notation 

k = 2Tr{m,n), x = {u,v)~^, dx = du dv, = {—iii,v,iii,u)~^ , F = {Fp, F^)~^ (B.3) 

The matrix elements with sine are defined accordingly. 


The elements of matrix N 22 in (17) 


N22 — 


( hii 

di2 

hic^ 

/il^/ 

\ 

^21 

^22 

h2c' 

~h2s' 


del 

hc2 

hcA 

hes' 



— hs2 


hss' 

J 


are dehned as 


dec' (^k^ k ^ 

hiAF) 

h^Ak') 


1 f • {k'R2A 

47r2 J I r2,^ X r2,„ | 

S 2 

1 f ^ r2,^ • {k' R2,x') i,, ^ 

— / dx— -— cos(/cx) 

27r J \ Y2,v X r2,n I 

S 2 

1 f - ll2,u ■ {k' R2,x'^ 


cos(A;'a;)) 


hii — 


— dx 

2vr J I T 2 ,v X r 2 ,„ , 

S2 

f , ll2,v ■ ll2,v , /■ , ll2,v ■ ll2,% 

I dx- -r, hi 2 — / dx-. 


S2 


X r2,t. 


S2 


h2i — 


dx-, - — - ^22 = / dx- 


S2 


r2,v X r2,^. 
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r2,v X r2,, 

ll2,u • r2,i, 
Ii2,v X r 2 ,, 


(B.4) 


(B.5) 



















and the matrix elements with sine are defined accordingly. 
The matrix Mi^ (17) is defined by[^ 


m,5 = 


/ 

0 

0 

0 

-F'p 

0 

0 

0 

\ 


0 

0 

0 

-F!p 

0 

0 

0 



0 

0 

0 
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0 
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0 

0 

0 

0 

0 

-\kF5k,k' 
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0 

0 

0 

0 

0 

0 





0 

0 

0 

0 

0 

0 



0 

\kF5k,k' 

0 

0 

0 

0 
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V 

0 

0 


0 

0 

0 

0 

/ 


(B.6) 


The matrices i = 1, 2, 3 in (26) are obtained by discretizing the vacnnm energy term 

(B.7) 


u 

'-'si 

h' 

'^s2 

U 

'^ss' 

lA 

"sc' 

Y^F'p (^1, 

i - l)i^r 

0 

0 

Yl 

hi 

"c2 

'^cs' 

F 

'^cd 


with 




If If 

/ dx sm{kx) d\{u,v), &s 2 (^) — sin(fca;) d\{u^v) (B.8) 

dvr ./c, dvr ja 


d\{u,v) = {Fp + Fp ri^t,) ■ -I- du'dv'riy x ^ — 27r Fp 


Si 


|ri - r;|3 


d\{u,v) = {Fp + Fp ri_^) ■ f du'dv'rip> x z = 2, 3 

Jsi Ti ~ 

dl{u, v) = {Fp + Fp ri_^) ■ -f du'dv'riy x + 2n Fp 


(B.9) 


5i 


|ri-r;|3 


d\{u,v) = {F'p + Fp ri_^) • [ du'dv'vi^u' x z = 2, 3 

Js, Ti ~ I'd 


and with F = {Fp, FpY 


hY{k,k') = —-{kF) / dx cos{kx) al,{u,v,k') 

2 Js, 


hY{k,k') = —-{kF) / dx sm{kx) ap{u,v,k') 

2 Js, 


(B.IO) 


h\p{k,k') = —-{kF) (n 6k,k' + J dx cos{kx) al,{u,v,F) 


bY{k,k') = -(fcF) I^TT 4,fc'5i,i + J dx sm{kx) al,{u,v,k') 


note: ^k,k' — ^m,m' 
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with 


f dx'{Yiy X riy) 

(n - r;; 

|ri - r;|' 

/ da:'(riy x ri,„/) 
JSi 

(r* - r[] 

ki - r'll' 


— sin(fc'a;') 


sin(fcV), i = 2,3 


and al,{u,v,k') accordingly. 

C The subtraction method 


(B.ll) 


The elements of the self-indnctance matrix Mu (B.l) and (B.2) are Fonrier integrals over 


the plasma-vacunm interface. The integrands are singnlar so that a standard nnmeri- 
cal Fourier transform is not possible. The problem is solved by applying a subtraction 
method: An analytically integrable function with the same singular behaviour is sub¬ 
tracted. The analytical integral is then added to the numerically Fourier-transformed 
regularized function. 


The subtraction method will be demonstrated by treating a term of (B.2) 


1111 

g{m,n;m' ,n') = J du J dv J j dv' g{u,v,u',v') e 

0 0 0 0 


• n' ^ 27 vi( 7 nu+nv—m'u'—n'v') 


with 


g{u,v;u\v') = 


r,, ■ r,,/ 


r(n, v) — r{u', v') 


(C.l) 


(C.2) 


Expanding g at the singularity 


/ X X X X 

r - r = Yu du + v^ dv + — + Yuv ou dv + Yuv — 

with du = u' — u, dv = v' — V and replacing du, dv by tan(7r du)/7i, tan(7r dv)/7i one gets 
a periodic function which can be Fourier-transformed analytically with respect to u' and 


gsing{u,v;u' - u,v' - v) (C.3) 

_ ^ _ 

-^/r^ tan^(7r(M'—u)) -|- 2 r„-r.y tan(7r(M'^)) tan(7r(n'^)) -|- tan^(7r(n'^)) 

With the analytically computable integrals (see appendix D) 
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1 1 

I{m,n;a,b,c) = du dv 


TT e 


2'n-i{mu-\-nv) 


0 0 

one finds for the singular integral 

1 1 


a tan^(7rM) + 2 6tan(7TO) tan(7rn) + c tan^(7rn) 


(C.4) 


gsingim, n; m', n') = J du J dv rl I (m', n'; r^, r„ ■ r^, r^) e“ 
0 0 

The regularized part can be Fourier transformed numerically 

greg{m,n]m\n) = 

1111 


(C.5) 


(C.6) 


= j duj dv j du' j dv' {g{u, v; u', v') - gsing{u, v; u' -U,v' - v)) V'-n'.') 

0 0 0 0 

so that the Fourier transform of g{u, v; u', v') is given by 


g{m, n; m', n') = greg{m, n; m', n') + gsing{m, n; m', n') 


(C.7) 


A second type of singular integrals appears in a part of the vacuum energy matrix 


(B.7) -(B.ll). The subtraction method will be demonstrated for a term in (B.ll) 

1 1 

h(n,n;m',n') = [ du' [ dv'h{u,v;u',v') 


0 0 


with 




I r{u, v) — r(u', v') ^ 
the expansion of the integrand at the singularity is given by 

hsing{u,v;u' -U,v' -v) = 


(C.8) 


(C.9) 


(C.IO) 


7 r(r^ X r„) tam( 7 r(n'— m)) + 2 r„^tan( 7 r(M'—n)) tan( 7 r(n'^)) + tam( 7 r(n'^)) 


2 (r^ tan^(7r(M'^)) +2 r„-r^ tan(7r(M'—n)) tan(7r(n'^)) + tan^(7r(n'^)))3/2 

With the analytically computable integrals (see apendix D) 

K{m, n; A, B, C, a, 6, c) = (C.ll) 

1 1 

f r ( 4 tan^lTT?/,'! -I- 2 R tan iTr?/,') tan iTrid d tFtn’^f'iruAp^'^drnu+nv) 

= 71 du dv- 


0 0 


(a tan^(7ra) + 2 6tan(7™) tan(7rn) + c tan^( 7 rn))^A 
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one gets for the singular integral 


hsing{u, V] m', n) = K{m', n; A, B, C, r^, ■ r^, r^) (C.12) 

with A = (r^ X r„) ■ 5 = (r^ x r„) ■ r^v, C = (r„ x r„) ■ 


The regularized integral can be Fourier transformed numerically. 

1 1 

hreg{u,v;m ,n) = j du J dv'{h{u,v;u',v') — hsingiu,v;u' — u,v' — v)) e 
0 0 

For the Fourier transform of h{u,v]u',v^) one gets 

h{u, v; m , n) = hreg{u, v; m', n') + hsing{u, v, m , n') 


2m{rn'u'+n'v') 

(C.13) 
(C.14) 


D The regularisation integrals 


In appendix C a subtraction method to regularize singular Fourier integrals has been 
presented using the following analytically computable integrals 
1 1 



Imn =71 dudv 


^27ri{mu-\-nv) 


0 0 


(atan^(7ra) + 26tan(7™) tan(7rn) + ctan^(7rn))2 


(D.l) 


and 


1 1 


, , (Atan^lTTu) + 2S tanf?™) tanlTTu) + Ctan^(7rn))e^”(™'“+”'^^ , 

Km,n = 77 du dv -3- 

(atan^(7ra) + 26tan(7rM) tan(7rn) + ctan^(7rn))2 



0 0 


with ac — b‘^ > 0. The integrals K^n can be obtained by deriving the integrals Imn with 
respect to a,b,c ■. 

o o O 

Kmn = -2(41— +B- + C-)Imn • (D.3) 

oa ob oc 

To compute the Imn a generating function X is introduced: 


X= ImnS^t^. 

m=0,n=0 


(D.4) 


Summing up the power series, one obtains X in closed form: 
1 1 

X = 7r 



du dv 


0 0 


(1 — se^™)(l — fe^™)(atan^(7™) + 26tan(7™) tan(7rn) + ctan^( 7 rn ))2 


(D,5) 
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With the variables (r, y) ■. {y = tariTTW, ry = tannv) one gets for X 
+ 00 +00 



I= — I I drdy 


— oo —oo 


y + ia y -ij \ry + i/3 ry - i J {^a + 2br + cr‘^)^ 


(D. 6 ) 


with 


1 -s ^ 1-t 

a = -, p =-. 

■’ ^ 1 + t 


1 + s' 

Integrating X with respect to y, one obtains 


I = -lir 


2 / * 


1 1 
+ 


/3 + ar 1 + r J — 2br + cr^) 2 
1 ^ 1 \ 1 


1 + ar p + r J (^a + 2br + cr^) 2 


and 


hiih) 


1 


Ij a{'y^/c + ac + (3b) 
7 I3{y^/a — ab — (da) 


(D.7) 


(D. 8 ) 


(D.9) 


with 7 = A/a+~+~ 26 c 7 ++c^ 

With the substitution x : a; = (1 — r)/(l + r) the function X can be written as a sum of 
four terms: 

X = t) + h’''(0, 0) + h“(s, 0) + h“(0, t ), (D.IO) 


with 


h^(s, t) 


(1 + s)(l + t) r dx 

2 y {1 — st — {s — t)x){a^+ 2dx + a^x'^)^ 


(D.ll) 


and 


— a 2b c, 

a~ = a — 2b + c, 

d = c — a. 


The Iran are then obtained by expanding X again as a power series in s and t. One starts 
by expanding the functions 


h^{s, t) 


1 (1 + s)(l + t) f S -t \ ^ , 

2 I-St ^ 


(D.12) 
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with 


+1 


Tr = dx 


X 


-1 


(a'f + 2dx + a^x ‘^)2 


Expanding the further, one hnally obtains 


{ C+n+ci-ln+ctin-l+C^-ln-l: 771 >l,7l>l 

cio + C:!n-10+ C;^0 + C^n-lO > 771 > 1, n = 0 

Con + Con-i+ c^n + Con -1 , form = 0,n > 1 

Coo+ Coo + Coo + Coo ,form = 0,n = 0 


with 


= 

^mn 


^ m+n+|m-n| 


E 

£=0 


+ ! T^_ri\+2i 


2 ( 


m+n—|m—n| 


i)\{\77l - 77\ + i)m 


The integrals can be calculated by using a recurrence relation 


(D.13) 


(D.14) 





y/c + c±b 

y/a — a^b’’ 



^(2(yc+ (-l)'ya) - (2f- l)(c- o)r,*i 


(D.15) 

{i — 1) ^£^ 2 ) for ^ > 2 . 


The Fourier integrals Kmn follow from the Imn by differentiation. One gets the same 
formulas (D.14) as for the Imn replacing the integrals by the integrals S^, which are 
given by 


c^± 



{A^ + 2D X + A^x‘^) 
(a^ + 2dx + 


(D.16) 


with 


A+ = A + 2B + C, 
A- = A-2B + C, 
D = C -A. 


As for the T'^, recurrence formulas can be derived for the 5*^, but the Sf can also be 


£ > 
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expressed in terms of the by appropriate partial integration. One obtains for S'/ 

(a + 26 + c)(ac — 6^)5'/ = 

+ 2B + C)(ac — h^) + l{ki{a + 26 + c) + ^ 2(0 — a))^ T/ 

+ ^{kl{c — a) + k 2 {a — sh + c)) (D.17) 

(c + 6 )fci + (c - 6)^2 . £ (a + 6 )/ci - (a - 6)^2 

with 

ki = (a + c)B — (A + C)b, 

k 2 = C{a + b) — B{c — a) — A{c + b). 

The formula for is obtained by replacing 6 by —6 and B by —B, and by 

T- T- 


The forward computation of the recurrence relation T/ gets unstable if 6 < 0, and for 
if 6 > 0 [32j. Considering the homogeneous equation 


£[T+| ^ T. 


+ 


with the ansatz TT 


{i + l)a 

{i + 2)a+ ' (/T 2 V 

one gets the characteristic equation 


T/ = 0 


{i + 2)a 

with the solution for large £ —>■ 00 


^2 (2£ + 3)d ^ _^{^+ l)a 


(7 + 2)a^ 


= 0 


d .\/a c — b^ , 

t± = —- ± t ---and 

a+ a+ 


U \= — 


(D.18) 


(D.19) 


(D.20) 


For 6 < 0 the moduli are | |> 1 so that the forward recursion is unstable. In that 

case the T/ can be obtained by solving the recurrence relation in the backward direction. 
First one generates a solution of the nonhomogeneous equation by using the equation in 
the backward direction with starting values 


zn{N + 1 ) — zn{N) — 0 . 

Then one generates a solution yN{(^) of the homogeneous equation C[y{i)] = 0 by using 
the equation in the backward direction with 

yN{N+ 1) = 0, yN{N) = l. 
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Then the solution wn{() is obtained by 


with 


WN{i) = KN)yN{tj + ZN{tj, ^<l<N + l 

- zn{Q) 


(D.21) 


\{N) = 


2/w(0) 


For N ^ oo the wn{() converges to . 


For large values of m,n an asymptotic expansion for the Fourier integral (D.l) 

^2'Ki{rau-\-nv) 


1 1 
2 2 


Inin ^ 



du dv 


1 _1 
■ 2 2 


(atan^(7rM) + 26tan(7rM) tan(7rn) + ctan^(7rn ))2 


can be derived. 

With the substitution 

and 


Tiu = r cos(p, TTV = r sin(p 


m = Xcos(po, n = Asin(po, A"'= 

27r R{<p) 

Imn = j dip j dr g{r, ip) 

0 0 


one gets 


with 


(D.22) 


(D.23) 


(D.24) 


fir.ip) 


2r cos((p — ipo) 

r 

7r(a tan^(r cos(p) + 26tan(r cos ip) tan(r sin(p) + ctan^(r sin(p)i/2) 


and 

R(cp) = J 2 |cos(c,)| for I _ lTr/2 |< vr/4 j ^ ^ 

[ 2|sinfo)| /-i,3 

The procedure used is based on the method given in [33]. In the present case the main 
contribution to the asymptotic expansion comes from the stationary point, where the 
derivatives ^f{r,ip) = 0, -^f[r,ip) = 0 vanish. The stationary point is found to be 


r 

cos((p - ipo) 


0 , 

0 ip - IPQ = -. 


(D.25) 
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Defining 99 ' : (^' = 99 — (^0 — f one gets with g{r, if') = g{r, ip) and f{r, ip') = /(r, ip) 


R 


Imn= / dr / dp' g{r,p') 


(D.26) 


-R -f 


with 


¥^') = 


7r(a tan^(r sin((^' + (^o)) — 26tan(r sin((^' + (^o)) tan(r cos((^' + (^o)) + ctan^(r cos((^' + (^o))^^^ 


There has been made use of the fact that g{—r, p) = —g(r, p) and g{r,p + tt) = g{r, p) . 
For R it is sufficient to choose R > 0 small but finite. 


Substituting p' hy t = sin((y9) one gets 


R 1 

dmn = J dr j dt G{r, t) e 
-R -1 


2i\rt 


(D.27) 


with 


G{rR) = 


7r(l — tan^(rTs) — 2b tan(rTs) tan(rTc) + c tan^(rTc))^/^ 


and 


Ts = t cos po + Vl — R sin(y9o, 

Tc = —t sin 990 + y/l — R cospo- 


(D.28) 


Substituting r = x — y, t = x + y one gets 


Imn = ‘2 dx dy G{x-y,x + y)) 


(D.29) 


The asymptotic expansion obtained for this case in [33] is given by 


Imn 


C,., = 


-2n £ 


c,, e 


u=Q 

\i+l 


(2A) 


i2+l 


ly ! ^^ 

i=o 

1 ■ 3 ■ ■ ■ (2£ - 1) 

229+1 \ J J {2^ + 2]){2^ + 2i-2)■■■{2] + 2) 


(- 1)^+1 / 2 j 


F-k = - — 


d' 

dx^ dx^ 


G{x-y,xRy) 


cc=0 

j /=0 


(D.30) 
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(D.31) 


The coefficients Cj^^-j can be written also in closed form 

^ (-l)^-.+i (2^-)! z ^ (2Z/-2J)! 

22i+i ]\]\N' 2"-i(z/-l)!’ 


N = 


2 ^-V! 


SO that one finds 


Un 


47r ^ (-1)^-^' 

j/l (^ 2 A)'^+^ 22 *^+! j\ {y — j)\dx‘^^ 


u\ 

-^G{x-y,x + y) (D.32) 


or 


Imn ^ ^ 


dvr 


^ii/7r/2_ 


z/! ( 2 A)^+^ 22^+1 ^5x2 5|/' 


^2 d^\", 

^ G(x-?/,x + 2/) 


or 


dmT7. 


and finally 


dTT 


^iv-wy _ 


^ z/! ( 2 A)‘^+i" 22^+1 

i /=0 ^ ' 


dr ^ dt 


■| + |) I «(’■-*) 


r’ = 0 
f =0 


^mn ^ ^ 


27r 


1-4) G(r,t) 


^ z/! ( 2 A)^+^ \dr dt 
The two leading terms of the asymptotic expansion are given by 

'G{r,t) G{r,t)rrtt 


r=0 

t=0 


(D.33) 


I =71 

J-mr) !' 


A 


8A3 


r=0 

t=0 


(D.34) 


Here it has been made use of the fact that G{r,t) is a function of r‘^, so that uneven 
derivatives vanish for r = 0 . 


Expanding G(r, f) at r = 0, f = 0 up to second order in r and t one gets 

G(r,i) = ^ 


7r^(l -t 2 )(Q,(t) + f3{t) r 2 )’ 
a(t) = a Tg — 2b TsTc +c T^, 

m = l(aTt-b{T.T! + T,Tl)+cTt). 


(D.35) 


Using (D.28) one gets for a{t) and /3(t) up to second order in t 

ait') = OiQ T Oil t T (X 2 1 /d(t) — bo bi t b2 


(D.36) 


where the coefficients Oj, b^, i = 0,1, 2 depend only on the m, n-harmonics (D.23) 

n 


m 


cos(po = -r, sin(po = w; ~ \/m2 + n2, 
A A 
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and are given by 


Qq = a sin^(ipo)—2b sin((^o) cos( 99 o) + c cos^((^o)) 

oi = (a — c) sin(2(^o) — 26 cos(2(^o), 

02 = {a — c) cos( 2 (y 9 o) + 26 sin(293o), 

2 

bo = -{a sin^i^po) - b sm{^po) cos{^po) + c cos^i^po)), (D.37) 

2 

61 = -(2 (a sin^((^o)-ccos^((^o))sin(2(^o)cos(2(^o)), 

4 2 

62 = (a + c) sin^(293o) - ^(® sin‘^(</9o) + c cos‘^(93o)) + sin( 2 (y 9 o). 

o o 

For the two leading terms of the asymptotic expansion one obtains 


_ 1 1 1 
f TD.rj, ^ ,- ~1~ 


8A3 I 


60 + 2 62 gQ- 2^0 + Q -161 ^ 15 0^69 

O E + 


\/ho 


4 


(D.38) 


The asymptotic expansion of Kmn (D.2) can be obtained from Imn using (D.3) 


K = - 

^ mn ^ ;_3 


1 ^0 3^0 / 60 + 2 62 0260 + 0161 350^60 


A y/oo 
1 

4A3 


+ 


8A3 


i/ho 


-5- 


i/ho 


+ 




(D.39) 


Bq + 2 i ?2 ^^260 + ®2 -Bo + ^161 + Oi-Bi ^ 15 2oiAi6o + (i\Bo 

o ^ E I : 


a/Oq 


-y/ho 


i/ho 


where the expressions for Ai,Bi,i = 0,1,2 are obtained from (D.37) by replacing a,b,c 
by A,B,C. 
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